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1  Introduction 


There  has  been  a  surge  of  interest  in  application  involving  gradiometer  data  recently,  particularly 
gradiometer  inverse  problems.  One  main  areas  of  application  surrounds  gradiometer  inverse  problems 
focused  on  the  European  Space  Agency's  geodetic  satellite  mission,  GOCE,  and  the  gravitational 
gradient  observations  produced  along  its  orbits.  This  inverse  problem  is  interested  in  producing 
spherical  harmonic  series  representations  of  the  Earth's  gravitational  potential  (Freeden  and  Nutz  2011; 
Murbock  et  al.  2011;  Novak  and  Tenzer  2013).  The  second  main  focus  area  is  on  geophysical 
prospecting  problems  now  feasible  due  to  recent  improvements  in  airborne  and  land  gradiometers,  such 
as  those  available  from  Lockheed  Martin  (Difrancesco  2007).  Additionally,  emerging  systems  based  on 
atom  interferometry  show  promise  at  increasing  instrument  accuracy  by  an  order  of  magnitude 
(McGuirk  et  al.  2002).  Also  advances  in  algorithms  aimed  at  modeling  the  gradients  from  local  terrain 
and  improving  the  likelihood  of  solving  the  local  prospecting  inverse  problem  are  presented  in  Jekeli 
(2012)  and  Uzun  and  Jekeli  (2015). 

The  inverse  source  problem  for  the  gradiometer  tensor  can  be  stated  generally  as  follows:  given  a 
gradiometer  tensor  field,  extract  information  about  the  unknown  object  from  which  it  was  generated.  In 
practice,  information  about  the  unknown  object  is  determined  by  identifying  model  parameters  that 
generate  gradiometer  terms  at  the  surveyed  locations  that  are  a  close  fit  to  the  observations.  This  may  be 
done  by  some  estimation  process  like  least  squares  or  by  trying  many  models  (forward  modeling)  and 
keeping  a  rbest  fittingO  one.  This  model  is  then  assumed  to  have  something  in  common  with  the 
unknown  object  I  location,  mass,  etc.  In  Anderson  (2011),  it  was  shown  that  the  gradiometer  inverse 
problem  reduces  to  the  inverse  problem  of  the  potential.  Mass  distributions  of  different  mass,  size,  and 
location  can  produce  very  similar  external  potential  fields.  Given  this  nature  of  the  inverse  problem, 
combined  with  instrument  and  environmental  noise  found  in  practical  applications,  the  inverse  solution 
is  often  wrong  in  fundamental  ways.  That  is,  the  inverse  problem  in  the  typical  geophysical  prospecting 
setting  is  usually  considered  unstable  unless  certain  parameters  are  highly  constrained.  In  this  paper, 
we  present  a  new  theory  that  specifies  when  an  inverse  solution  model  would  have  the  same  center  of 
mass  and  total  mass  as  the  unknown.  The  total  mass  estimate  includes  estimating  regions  with  both 
negative  and  positive  density  contrast.  In  practical  problems,  this  can  help  mitigate  the  pervasive 
instability  problem  that  can  produce  false  positives/negatives  typically  caused  by  not  correctly 
determining  the  unknown's  true  center  of  mass  location  and  total  mass.  In  addition,  we  derive  a  bound 
for  the  center  of  mass  location  and  the  total  mass  for  strictly  positive  (or  negative)  bounded  mass 
anomalies  as  a  function  of  the  supremum  norm  of  the  differential  curvature. 

2  Mathematical  Basis 


Let  ju(y)  be  a  mass  distribution.  The  potential  function,  U,  generated  by  jd  is  defined  by  the 
Lebesgue  integral 


|x  -  y| 


l 


where  x  is  a  point  in  space,  y  is  a  point  in  the  mass  distribution,  and  G  is  Newton's  gravitational 
constant. 

The  Hessian  of  the  potential,  H(U ),  is  the  tensor  field  whose  components  in  rectangular  coordinates 
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are  Utj  = - .  It  is  often  convenient  to  express  H in  matrix  notation. 
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where  V  = 


'  d_  d_  _S_"| 
v  dx  ’  dy  ’  dz  J 


is  the  gradient  operator. 


The  matrix  in  Eq.  2  is  symmetric  and  satisfies  LaplaceQ  equation.  Thus,  H  has  only  five 
independent  terms  and  is  a  symmetric  order  2  tensor  field  outside  the  support  of  the  mass. 

The  Partial  Tensor  gradiometer  produces  observables  comprised  of  a  combination  of  Hessian  terms 
restricted  to  the  horizontal  plane,  referred  to  as  the  in-line  and  cross  gradients  and  defined  as  P=  Uxx  - 
Uyy  and  Q=2Uxy  respectively.  There  is  also  a  Full  Tensor  gradiometer  that  supplies  H  as  the  observable, 
but  our  work  focuses  on  the  Partial  Tensor  system. 
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The  Gradiometer  Tensor  generated  by  the  mass  anomaly  ju,  denoted  GTM,  is  the  second-order  tensor 
field  in  the  xy-plane  defined  by 

GT M  =  P(dx  <8>  dx  —  dy®  dy )  +  Q(dx  ®  dy  +  dy  <8>  dx).  3 

Note,  the  gradiometer  tensor  is  a  symmetric  and  trace  free  (trace  =  0)  co variant  tensor  field  of  order  2 
on  the  xy-plane.  Thus,  Eq.  3  can  be  written  in  matrix  notation  as 

GTM ={P  4 

U2  ~P) 

From  the  gradiometer  tensor,  the  Differential  Curvature,  which  is  the  scalar  field  DC  in  the  plane,  is 
readily  calculated  as 


DC  =  JpT7q2  5 

The  points  where  the  differential  curvature  is  zero  (i.e.,  P=Q=0)  are  exactly  the  singularity  points  of 
the  gradiometer  tensor.  It  was  shown  in  Anderson  (201 1)  that  the  value  of  DC  does  not  change  with 
coordinate  system  rotation  or  translation,  which  makes  the  differential  curvature  a  convenient 
observable  in  practical  applications. 

A  simple  calculation  shows  that  the  eigenvalues  of  the  gradiometer  tensor  at  a  given  point  are  ±DC. 
Since  the  gradiometer  tensor  is  symmetric  at  points  where  it  is  non-singular,  the  eigenvectors  associated 
with  +DC  and  V  DC  are  mutually  perpendicular. 

Define  the  Line  Field  A  associated  with  GT1  as  the  line  field  of  GT p  such  that  A(x,y)  is  the 
eigenspace  of  GP1  at  the  non-singular  point  (x,y)  corresponding  to  the  positive  eigenvalue  DC(x,y).  It  is 
convenient  to  describe  the  line  field  by  giving  the  angle  odx,y)  that  A(x,y)  makes  with  the  x-axis.  This 
angle  is  called  the  azimuthal  angle  of  GT1  at  (x,y).  The  formula  for  the  azimuthal  angle  of  a 
gradiometer  tensor  in  terms  of  its  components  is  the  following: 

0.5  tan -1  —  P*  0 

a  =  ^  6 

+  —  P  =  0 

l  2 

The  arctangent  in  Eq.  6  is  computed  such  that  the  signs  of  P  and  Q  are  accounted  for,  resulting  in  an 
angle  that  ranges  from  - '  to  + Thus,  a  is  on  the  range  of  [-  72,  +  72]. 

The  gradiometer  tensor  at  any  point  is  completely  determined  by  knowing  its  eigenspace.  Thus,  by 
knowing  DC  and  a,  we  know  the  tensor  itself.  The  reason  for  using  DC  and  a  is  that  it  allows  one  to 
better  visualize  the  gradiometer  tensor  field. 

Now  we  apply  index  theory  results  to  the  gradiometer  tensor  line  field,  and  show  that  the  global  index 
value  of  the  line  field  can  provide  useful  information  about  the  unknown.  Background  on  index  theory 
can  be  found  in  section  6.8  in  Strogatz  (1994)  and  is  summarized  below. 

If  C  is  a  closed  curve  in  the  plane  which  does  not  pass  through  any  singularities,  then  the  index  of  the 
curve,  Ic,  with  respect  to  the  field  is  defined  by  the  line  integral 

I  r  =  —  [da  1 

C  2n  J 
c 

Geometrically,  the  index  measures  the  net  rotation  of  the  line  field  as  one  goes  around  the  closed 
curve.  The  index  of  an  isolated  singularity  is  defined  to  be  the  index  of  any  simple  closed  curve  with 
the  usual  orientation  that  surrounds  only  that  singularity.  Some  properties  of  the  index  follow. 

1)  For  any  closed  curve  C,  Ic  is  an  integer  multiple  of  1/2  which  is  positive  if  the  net  rotation  is 
counterclockwise  and  negative  if  clockwise. 

2)  For  a  simple  closed  curve,  C,  with  the  usual  orientation  that  surrounds  only  isolated  singularities,  Ic 
is  the  sum  of  the  indices  of  the  singularities. 

The  line  field  index  is  said  to  be  a  global  index,  /G,  if  the  curve  C  surrounds  all  of  the  singularities  of 
the  field. 

Fig.  1  shows  the  line  field  for  three  point  masses  at  (-2,  -1,  -2),  (0,  -2,  -2),  and  ( 1, 1,-2)  with  masses 
0.5,  -1.5,  and  2  respectively.  Areas  in  which  P<0  are  indicated  with  a  blue  'x'  and  areas  in  which  Q<0 
are  indicated  with  a  red  'o'.  Singularities  occur  at  the  intersections  of  the  P=  0  and  Q= 0  contours.  Each 
cycle  through  (P>0,  Q>0),  (P<0,  Q>0),  ( P<0,Q<0 ),  (. P>0,Q<0 )  represents  a  180  degree  rotation  of  the 
line  field  and  an  index  of  0.5.  Progressing  through  this  sequence  in  the  opposite  direction  indicates  a 
negative  singularity  with  an  index  of  -0.5  Using  the  perimeter  of  the  plot  as  the  curve  C  produces  two 
complete  cycles  and  an  index  of  1,  as  predicted  from  Theorem  1  for  this  unbalanced  system.  Note  that, 
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if  curve  C  were  a  circle  with  radius  10  centered  at  the  origin,  it  would  exclude  the  negative  singularity 
at  (0,-15)  and  indicate  an  index  of  1.5.  This  result  illustrates  the  importance  of  ensuring  C  encloses  all 
singularities  when  determining  the  global  index. 
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Fig.  1  Example  line  field  for  a  system  of  3  point  masses 

A  mass  anomaly  is  the  difference  between  the  unknown  object  and  an  assumed  model  (the  inverse 
solution).  For  a  given  mass  anomaly,  the  places  not  modelled  correctly  can  be  considered  regions  of 
positive  or  negative  mass.  Given  the  mass  anomaly,  //,  define 


m 


as  the  positive  and  negative  parts  of  the  mass  anomaly  respectively.  Define  the  centers  of  mass  with 
respect  to  the  x-axis  for  the  positive  and  negative  parts  of  the  mass  anomaly  as 


-+  _  f  _  f  %dju 

m+  ^  m~ 
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The  centers  of  mass  with  respect  to  the  j-axis  for  the  positive  and  negative  parts  are  similarly  defined. 


Theorem  1:  (1)  If  m +  ^  m  then  the  global  index,  IG,  equals  1.  (2)  If  m+  =  m  and 

(x  ,y  )  *  (x  ,y  ),  then/(7=3/2. 

Corollary  1.1:  If  IG  is  not  equal  to  one,  then  m+  =  m\  and  if  IG  is  greater  than  or  equal  to  2,  then 

(x+ ,  y+ )  —  (x  ,y  ). 


Theorem  1  and  its  corollary  from  Anderson  (201 1)  state  that  if  the  gradients  from  the  mass  anomaly, 
(i.e.,  the  residual  gradients  calculated  from  the  difference  of  the  gradients  generated  by  the  inverse 
solution  model  and  the  gradients  produced  by  the  unknown  object)  have  a  global  index  02,  then  the 
inverse  solution  has  the  same  total  mass  and  center  of  mass  in  (x,y)  as  the  unknown.  The  conditions  of 
Theorem  1  are  illustrated  in  Fig.  2.  These  plots  show  the  line  fields  produced  by  a  prism  and  a  cylinder. 
Formulas  used  to  calculate  the  gradients  from  a  prism  or  cylinder  can  be  found  in  Dransfield  (1994)  and 
Romaides  et  al.  (2001)  respectively.  In  Fig.  2a,  the  cylinder  mass  is  1.3  times  the  prism  mass  and  it  is 
offset  3m  along  the  x  axis,  resulting  in  an  index  of  1.  In  Fig.  2b,  the  masses  are  the  same,  but  the 
cylinder  is  still  offset  by  3m,  resulting  in  an  index  of  1.5.  Finally,  in  Fig.  2c,  both  the  masses  and 
centers  of  mass  of  the  two  objects  are  the  same,  resulting  in  an  index  of  2.  The  far  left  singularity  in 
Fig.  2a  has  an  index  of  -1/2,  while  all  other  singularities  have  an  index  of  +1/2.  This  produces  an  IG=  1 
in  accordance  with  Theorem  1. 
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Fig.  2  line  fields  for  a  cylinder  and  prism  with  (a)  different  mass  and  center  of  mass,  (b)  same 
mass  and  different  center  of  mass,  and  (c)  same  mass  and  center  of  mass.  The  circles  represent 
singularity  points  where  P=Q=0. 

Theorem  1  and  Corollary  1.1  also  apply  for  the  vertical  planes,  P(x,z)  and  P(y,z).  So  for  certain 
applications  where  the  anomaly  is  bounded  by  areas  below  both  a  horizontal  and  a  vertical  plane  then 
we  get  the  result  for  3D  center  of  mass. 

Corollary  1.2:  If  IG  >  3/2  in  any  two  of  the  three  planes  P(x,y),  P(x,z),  P(y,z)  then 
( x  ,y  ,  z  )  =  (x  ,y  ,  z  ),  and  m+  =  m. 

The  following  proposition  from  Anderson  (2011)  applied  with  Theorem  1  is  used  as  rationale  behind 
the  new  algorithm  introduced  in  the  next  section.  For  a  positive  integer  TV,  define  the  2N-gon  mass 
anomaly  by  placing  a  mass  with  sign  (-l)k  at  the  point  locations  (cos(k7i/TV),  sin(k7i/TV),  -1);  (k  =  0  e  2N 
-  1).  This  anomaly  has  positive  and  negative  masses  that  balance  and  have  the  same  center  of  mass. 

Proposition  1:  The  global  index  of  the  2N-gon  mass  anomaly  is  1  +  TV/2. 

The  line  field  for  a  2N-gon  with  TV=3  is  shown  in  Fig.  3  and  has  an  index  of  2.5. 


-15  -10  -5  0  5  10  15 

X 

Fig.  3  Line  field  of  a  2N-gon  for  N=3 

3  Results 

3.1  Minimum  Mass  until  Bifurcation  (MMB)  Algorithm 

Theorem  1  and  its  corollaries  are  extremely  sensitive.  Very  small  differences  in  the  mass  magnitudes 
or  center  of  mass  locations  result  in  IG< 2.  In  applied  settings,  we  assume  the  unknown  mass  is  large 
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relative  to  the  noise.  Additionally,  in  the  applied  setting  we  do  not  have  the  luxury  of  exact 
observations  or  surveying  great  distances  from  the  anomaly  to  get  a  global  index.  However,  even  in 
low  noise  settings,  this  noise  is  sufficiently  large  that  we  expect  a  local  index  of  1. 

The  pattern  produced  by  the  2N-gon  is  very  unlikely  to  occur  naturally.  This  makes  it  useful  in 
conjunction  with  Theorem  1  for  evaluating  how  accurately  a  model  represents  an  unknown  in  terms  of 
its  total  mass  and  center  of  mass.  We  place  a  model  prism  at  the  hypothesized  center  of  mass  of  the 
unknown  mass  and  place  a  2N-gon  just  below  the  measurement  plane  and  above  the  center  of  mass  of 
the  model.  We  then  increase  the  mass  on  the  vertices  of  the  2N-gon  until  its  distinctive  line  field 
pattern  emerges  at  a  selected  distance  from  the  center  of  the  2N-gon.  We  repeat  this  process  while 
varying  the  position  and  dimensions  of  the  model.  The  model  position  and  dimensions  that  require  the 
minimum  mass  on  the  2N-gon  to  produce  its  distinctive  pattern  represent  the  best  fit  to  the  unknown. 
We  call  this  smallest  2N-gon  mass  required  to  produce  the  pattern  the  minimum  mass  until 
bifurcation  (MMB). 

An  example  result  is  shown  in  Fig.  4.  The  unknown  is  a  rectangular  prism  1.25m  wide,  1.72m  tall, 
and  15m  long  with  a  density  of  2300kg/m3.  A  model  prism  of  the  same  length  and  density  is  placed  at 
the  unknown's  center  of  mass  and  its  width  and  height  are  varied  in  2cm  increments.  The  resulting 
surface  represents  the  MMB  required  on  the  2N-gon  for  each  combination  of  model  width  and  height. 
As  can  be  seen  from  the  plot,  the  overall  MMB  occurs  when  the  model  dimensions  match  those  of  the 
unknown.  The  plot  also  reveals  an  arc  of  local  minima.  These  points  represent  combinations  of  width 
and  height  which  produce  models  that  are  similar  in  total  mass. 


Width  (m)  Height  (m) 


Fig.  4  MMB  for  a  model  prism  that  is  a  near  perfect  match  in  dimensions  to  the  unknown.  The 
large  spike  occurs  at  the  correct  prism  dimensions. 

If  no  model  perfectly  matches  the  unknown,  the  MMB  may  occur  anywhere  along  the  arc.  For 
example,  if  the  closest  model  prism  width,  length,  and  height  have  a  1cm  error,  result  is  illustrated  in 
Fig.  5.  The  largest  spike  occurs  at  width  1.32m  and  height  1.63m. 


Width  (m)  Height  (m) 

Fig.  5  MMB  when  no  model  is  a  perfect  fit.  The  red  dot  represents  the  correct  prism 
dimensions. 

The  MMB  approach  can  produce  inverse  solutions  very  different  than  those  based  on  an  inverse 
model  solution  found  by  minimizing  the  difference  between  model  generated  gradients  and  the 
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observed  or  measured  gradients.  To  demonstrate  this,  in  the  next  example,  the  model  with  the 
combination  of  mass  and  location  that  minimizes  aDC ,  the  standard  deviation  of  the  difference  between 
the  DC  from  the  point  mass  model  and  the  observed  DC  of  the  unknown  prism,  is  compared  to  MMB 
point  mass  solution  to  the  same  unknown.  This  example  illustrates  the  false  positive  issue  often  caused 
when  the  model  has  very  different  shape  characteristics  than  the  unknown.  In  this  case  the  model  set 
can  be  considered  from  a  family  of  homogenous  balls  (point  masses),  while  the  unknown  is  a  prism. 

In  the  following  figures,  we  compare  inverse  solutions  for  the  unknown  prism  with  a  point  mass  and 
compare  the  results  of  the  MMB  algorithm  and  the  best  fit  approach  computed  over  increasing  distances 
in  the  plane  above  the  prism  center.  We  place  the  point  mass  at  one  possible  location  for  the  prism 
center  of  mass  and  compute  gdc ,  the  standard  deviation  of  the  difference  between  the  DC  from  the  point 
mass  model  and  the  observed  DC  of  the  unknown  prism.  We  then  place  a  2N-gon  above  the  point  mass 
and  determine  the  MMB.  This  process  is  repeated  for  all  combinations  of  mass  and  position  and 
identify  the  combinations  that  minimize  either  oDC  or  MMB.  For  both  approaches,  we  evaluate  the 
result  over  a  series  of  concentric  squares  to  examine  how  distance  from  the  prism  affects  the  results.  If 
the  two  methods  were  equivalent,  they  would  produce  the  same  result.  Both  methods  accurately 
determine  the  horizontal  center  of  mass  of  the  prism,  but  produce  different  results  for  its  mass  and 
depth.  For  a  1.0x20x1. 5m  prism  at  a  depth  of  12m  with  density  2300kg/m3,  Fig.  6  shows  oDC  as  a 
function  of  the  point  mass/prism  mass  ratio  and  the  horizontal  range  from  the  prism  center  of  mass  to 
the  edge  of  the  square  area  over  which  oDC  is  computed.  The  minimum  oDC  at  each  range  occurs  at  the 
largest  mass  ratio  tested.  Thus,  the  inverse  solution  that  most  closely  fits  the  observed  data  yields  an 
incorrect  mass  of  the  unknown  at  all  depths  tested. 


Min  aDC  (s'2) 


Fig.  6  Best  fitting  solution  using  minimum  oDC.  The  minimum  oDc  at  each  range  is  at  the 
maximum  mass  ratio  tested  and  decreases  as  range  increases. 

Fig.  7  shows  the  equivalent  result  when  using  the  MMB  approach.  As  the  range  increases,  the  point 
mass  that  minimizes  the  MMB  and  the  depth  error  remains  consistent. 


Fig.  7  MMB  as  a  function  of  range  and  mass  ratio.  MMB  at  each  range  occurs  when  model 
mass  equals  prism  mass 

Using  this  approach  of  evaluating  the  MMB  computed  along  a  series  of  increasingly  larger  concentric 
squares,  it  is  possible  to  determine  if  the  model  mass  is  too  large  or  too  small.  Based  on  these 
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observations,  we  conclude  that  there  is  greater  confidence  in  an  MMB  solution  that  is  stable  over 
increasingly  larger  ranges. 

The  magnitude  of  the  MMB  can  provide  some  insight  into  how  well  the  model  fits  the  unknown.  In 
Fig.  8,  we  replace  the  point  mass  used  to  generate  Fig.  7  with  a  cylinder,  which  is  a  much  closer  fit  to 
the  prism.  The  cylinder's  length  and  density  is  set  to  be  the  same  as  the  prism  and  its  radius  is  varied. 
Note  that  the  MMB  is  lower  at  the  correct  ratio  with  the  cylinder  than  it  was  for  the  point  mass. 

MMB  C o nt o u rs- Cylinder,  Model  Depth— 12 
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Fig.  8  MMB  as  a  function  of  range  and  mass  ratio  for  a  cylinder  modelling  a  prism 

In  a  theoretical  environment  with  perfect  observations,  the  MMB  solution  allows  us  to  find  the  correct 
mass  and  horizontal  center  of  mass  of  the  unknown.  However,  for  practical  problems  it  is  very  difficult 
to  separate  mass  and  depth.  Also,  as  was  illustrated  in  Fig.  5,  multiple  inverse  solutions  of  the  same 
mass  can  be  viable  candidates  for  the  best  fit. 

3.2  Mass  and  Depth  Bounds  from  DC  Values 

In  this  section,  we  derive  some  fairly  sharp  bounds  on  total  mass  and  the  center  of  mass  location 
(depth)  of  a  strictly  positive  (or  negative)  unknown  mass  that  are  a  function  of  the  greatest  magnitude  of 
the  observed  DC  values  (its  supremum  norm).  Our  approach  takes  into  account  gradient  measurement 
error  and  mass  modeling  errors,  assuming  the  unknown  is  contained  within  a  bounding  ball  and  is 
strictly  positive  or  negative.  If  the  unknown's  location  cannot  be  bounded,  then  neither  the  center  of 
mass  nor  the  total  mass  can  be  bounded. 

Set  a  unit  point  mass  at  depth,  D ,  on  the  z-axis  below  the  (x,  y)  computational  plane.  Then, 
symbolically, 


DC2  =(UXX-Uyy)2  +(2UX},)2 


9(x2  +  y2)2 
(x2  +  y2  +  D2)5  ' 
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Now,  find  the  maximum  DC  value  in  the  (x,  y)  plane.  Since  the  DC  surface  plot  is  rotationally 
symmetric,  we  can  further  simplify  Eq.  10  by  setting  y=0.  Then,  using  methods  of  Calculus,  we  take 
the  derivative  with  respect  to  x,  set  the  result  equal  to  zero,  and  solve  for  the  roots. 


SDCy=0  _3x(-3x2  +2 £>2) 

dx  1 

(x2  +D2)2 
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The  roots  of  Eq.  11  are  x  =  ±Z>V 2/3  .  Applying  these  roots  to  Eq.  10  with  y= 0  and  inserting  the 
gravitational  constant,  G ,  and  mass,  m,  we  get  a  relationship  between  the  supremum  norm  for  DC  and 
the  point  massa  depth,  D,  and  the  mass,  m, 

IMsup  12 

As  an  example  of  how  this  equation  would  be  applied,  consider  the  prism  used  for  Fig.  7.  The 
\\DC\\sup  for  this  object  is  1.012xl0'V2  and  its  mass  is  6.9xl04kg.  Using  these  values  in  Eq.  12  and 
solving  for  depth  produces  a  result  of  13.28m,  a  10.6%  error  for  an  object  quite  different  from  a  point 
mass. 


For  a  quick  and  simple  estimate  of  the  effect  of  observational  error  on  mass  and  depth  error,  we  take 
the  derivative  of  Eq.  12  with  respect  to  mass  and  depth: 


^IMsup 


18V3G  ,  5Ap,Gm 


dD. 
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To  demonstrate,  assume  the  observational  error  for  DC  for  the  preceding  example  is  +5xlO'10s'2. 
Using  the  mass  of  6.9xl04kg  and  the  computed  depth  of  13.28m,  solving  Eq.  13  for  dm  and  dD 
indicates  the  mass  may  vary  by  ±3. 14xl04kg  and  the  depth  may  vary  by  ±2.02m. 


Claim  1.  For  any  mass  distribution  generating  an  admissible  potential  in  the  ( x ,  y)  plane,  p,  outside 
ball,  B\  its  ILDCIIsup  is  bounded  below  by  ILDCIIsup  generated  of  a  uniform  ring  of  unit  mass  in  the 
equatorial  plane  of  B  parallel  to  p  and  bounded  above  by  the  ILDCIIsup  of  two  point  masses,  each  with 
mass  Vi  located  at  opposite  ends  of  the  diameter  of  B  (oriented  between  0  and  12  degrees  off  the  z-axis 
as  D  increases  from  1  to  infinity),  I IT)Glsup?ringo I LDCIIsup6 1 LDCIIsup ^  pts. 


Fig.  9  illustrates  the  results  described  in  Claim  1.  As  defined  in  Green  (1952),  an  admissible  potential 
is  generated  by  a  mass  distribution  within  a  unit  ball  where  all  the  mass  is  contained  in  the  ball,  the  total 
mass  is  1,  and  the  center  of  mass  is  at  the  center  of  the  ball.  We  solve  Eq.  12  for  D  and  compute  the 
depth  error  for  various  mass  distributions  using  their  calculated  DC.  The  horizontal  and  vertical  prisms 
used  for  this  plot  had  a  0.04m2  square  cross  section  and  a  length  that  placed  their  corners  on  the  surface 
of  the  ball.  The  errors  converge  toward  0  as  depth  increases  and  the  mass  distribution  more  resembles  a 
point  mass  from  the  plane. 


Depth  Estimation  Errors,  Sphere  Radius=1 ,  Mass=1 


Fig.  9  Depth  estimation  errors  of  various  distributions  of  a  unit  mass  in  a  unit  ball 
Heuristically,  we  argue  for  the  validity  of  the  Claim  1  by  citing  Theorems  2  and  3  from  Green  (1952). 


Theorem  2:  The  greatest  distance  from  the  center  of  mass  to  an  equipotential  surface  occurs 
when  the  generating  mass  of  the  admissible  potential  is  from  two  equal  point  masses  located  on 
the  boundary  of  the  unit  ball  where  the  two  points  are  aligned  with  the  furthest  point  on  the 
surface. 

Theorem  3:  The  closest  point  on  an  equipotential  surface  to  the  center  of  mass  occurs  when  the 
generating  mass  is  contained  in  an  equatorial  great  circle  with  the  closest  point  on  the  surface  on 
a  line  through  itself  and  the  center  of  mass. 


Empirically,  we  see  in  Fig.  9  how  the  ring  and  two  vertical  point  mass  distributions  bound  a  variety  of 
arbitrary  mass  distributions  contained  in  the  unit  ball.  Thus,  even  if  a  counter  example  is  found  and 
Claim  1  is  proven  false,  it  still  has  practical  utility  for  our  purposes. 

One  can  also  define  DC  in  any  Tangent  Plane  (TP)  to  the  equipotential  surface  normal  to  the  Z  axis, 
as  the  product  of  Uz  and  the  difference  between  maximum  and  minimum  curvatures  in  TP  (k1  and  k2 
respectively)  (Slotnick  1932). 


DC  =  Uz(kl-k2).  14 

We  also  note  that  the  ring  is  the  least  concentrated  mass  arrangement  on  the  equatorial  great  circle 
producing  relatively  small  Uz  values.  The  equations  for  the  gradients  of  the  ring  are  quite  complicated 
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(Lass  and  Blitzer  1983),  but  we  have  found  that  a  circular  array  of  point  masses  provides  a  good 
approximation.  Note,  the  maximum  DC  of  the  ring  and  two  vertical  points  configurations  occur 
respectively  farthest  from  and  closest  to  the  origin  of  the  computational  plane  as  illustrated  in  Fig.  10. 
These  facts  motivated  the  development  of  Claim  1. 


Fig.  10  DC  for  ring,  1  point,  and  2  point  masses  rotated  10  deg  from  the  vertical  at  6m  depth 

Fig.  11  shows  the  orientation  of  two  point  masses  on  the  boundary  that  produce  the  largest  ILDCIIsup 
values  as  a  function  of  depth.  The  orientation  approaches  vertical  as  depth  approaches  1  and  appears  to 
be  asymptotic  to  approximately  12  degrees  as  depth  approaches  infinity.  In  most  practical  applications, 
a  vertical  orientation  provides  an  effective  upper  bound. 


Angle  of  rotation  that  produces  max  DC  for  each  depth 


Fig.  11  Angle  from  vertical  at  which  two  points  produce  the  max  IIDCIIsup  as  a  function  of  depth 

We  now  consider  an  example  in  which  we  compute  bounds  for  the  unknown®  mass  and  depth  (center 
of  mass)  using  Eq.  12  and  Claim  1.  For  this  example,  we  use  a  prism  3x10x3m  at  a  depth  of  15m. 
Given  an  observed  maximum  DC  of  2.16xl0'9s'2,  we  assume  the  object  that  generated  this  DC  can  be 
enclosed  in  a  ball  with  radius  5.5m.  We  position  point  masses  at  the  top  and  bottom  of  the  ball  and  use 
Eq.  12  to  determine  the  mass  and  depth  combinations  that  produce  the  same  maximum  DC.  We  repeat 
this  process  using  a  ring  at  the  ball's  equator.  The  result  is  a  pair  of  curves  that  bound  the  depth  and 
mass  of  the  unknown,  represented  by  the  dashed  lines  in  Fig.  12. 

Reading  from  Fig.  12,  we  can  see  the  depth  range  attributed  to  the  inverse  solution  with  a  mass 
determined  from  Fig.  8.  For  an  object  with  an  MMB  solution  mass  of  2.07xl05kg,  the  depth  to  the 
center  of  mass,  D,  falls  in  the  range  -13.7m>D>-18.1m  based  on  the  mass  distributions  as  discussed  in 
Claim  1. 
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If  we  assume  a  maximum  error  signal  of  ±7xlO'10s'2  and  apply  this  to  the  DC  value  used  to  compute 
the  depth,  we  get  the  bounds  represented  by  the  solid  curves  in  Figure  12. 


x  10* 


Depth  (m) 


Fig.  12  Relationship  between  depth  and  mass  for  two  vertical  points  (lower  right)  and  a  ring 
(upper  left)  on  a  ball  with  radius  5.5m.  Dashed  lines  represent  bounds  with  no  measurement  error. 
Solid  lines  represent  bounds  with  7xlO'10s'2  measurement  error. 

The  depth  bounds  calculated  above  will  be  wider  for  shallower,  less  spherical  mass  distributions 
than  those  we  calculated  from  Eq.  12.  The  bounds  will  approach  each  other  as  depth  increases 
with  a  fixed  bounding  ball  as  all  distributions  will  then  resemble  a  point  mass  from  the 
observational  plane. 

4  Conclusion 

In  conclusion,  we  presented  an  algorithm  (the  MMB)  that  may  be  less  susceptible  to  false  positives  in 
geophysical  prospecting  applications  than  those  approaches  relying  solely  on  methods  selecting  inverse 
solution  models  that  most  closely  fit  the  observed  gradient  data.  This  is  true  since  the  MMB  algorithm 
produces  solutions  that  have  the  same  horizontal  center  of  mass  as  the  unknown  and  that  ensure  the 
anomaly  (difference  between  inverse  solution  mass  and  unknown  mass)  has  approximately  the  same 
amounts  of  positive  and  negative  mass.  Additionally,  we  presented  tools  for  bounding  possible 
locations  of  the  unknown  center  of  mass  and  total  mass  for  strictly  positive  (or  negative)  mass.  These 
bounds  alert  us  to  the  condition  when  the  size  of  the  unknown  mass  may  be  small  in  relation  to  the 
noise  in  the  particular  inverse  problem,  increasing  the  danger  of  an  ill-posed  problem  and  increasing  the 
likelihood  of  false  positives  or  negatives. 

While  this  analysis  assumed  a  dense,  regularly  spaced  array  of  observations,  future  work  will  examine 
the  application  of  the  approaches  presented  here  to  sparse  observations  sets. 
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